00100		SUBROUTINE LPC(AIFFY,SPT,NPTS,M,NSP)
00200	C	THE PARAMETERS IN THE LIST ARE DEFINED AS FOLLOWS:
00300	C	AIFFY←THE INPUT DATA ARRAY OF REAL NUMBERS...A TIME SERIES
00400	C	NPTS←THE NUMBER OF POINTS IN THE TIME SERIES
00500	C	M←   THE NUMBER OF FILTER COEFFICIENTS DESIRED.
00600	C	SPT← THE ARRAY IN WHICH THE REAL SPECTRUM IS RETURNED.
00700	C	NSP← THE NUMBER OF POINTS IN THE REAL SPECTRUM.
00800		DIMENSION AIFFY(1),CF(50),SPT(1),X(512),Y(512)
00900		TPI=2.*355./113.
01000	10	NP=NPTS-1
01100	C	CALCULATE THE DIFFERENCE SIGNAL
01200		DO 11 J=1,NP
01300	11	X(J)=AIFFY(J+1)-AIFFY(J)
01400	C	APPLY A HAMMING WINDOW TO THE TIME SERIES.
01500		DT=NP
01600		X(NPTS)=X(NPTS-1)*.08
01700			DO 12 J=1,NP
01800			T=J-1
01900	12		X(J)=X(J)*(.54-.46*COS(TPI*T/DT))
02000	C	CALCULATE AUTOCORRELATION COEFFICIENTS AND SOLVE
02100	C	FOR INVERSE FILTER COEFFICIENTS
02200	100	DO 105 J=1,NPTS
02300	105	Y(J)=X(J)
02400		CALL INVFL(Y,NPTS,M,RO,ERRN)
02500	C	Y ARRAY RETURNS 1.,A1,A2,...,AM,0.,....0.
02600		MM=M+1
02700		DO 125 J=1,MM
02800	125	CF(J)=Y(J)
02900	300	DO 303 J=1,512
03000		X(J)=0.
03100	303	Y(J)=0.
03200		MM=M+1
03300		DO 310 J=1,MM
03400	310	X(J)=CF(J)
03500	C	SHUFFLE DATA FOR CALC OF REAL TRANSFORM
03600		MM2=MM/2 +1
03700			DO 320 K=1,MM
03800			K1=2*K-1
03900			K2=2*K
04000			X(K)=X(K1)
04100	320		Y(K)=X(K2)
04200		NMOD=1
04300	323	IF (2**NMOD.GE.NSP) GO TO 326
04400		NMOD=NMOD+1
04500		GO TO 323
04600	326	IF ((2*MM2).LT.MM) MM2=MM2+1
04700		MOD=1
04800	330	IF (2**MOD.GE.MM2) GO TO 340
04900		MOD=MOD+1
05000		GO TO 330
05100	340	CALL FFTP(X,Y,NMOD,MOD)
05200		CALL RBITS(X,Y,NMOD)
05300		CALL USCRM(X,Y,NMOD)
05400	C	CALCULATE LOG OF RECIPROCAL INVERSE FILTER SPECTRUM
05500			DO 350 J=1,NSP
05600	350		X(J)=X(J)*X(J)+Y(J)*Y(J)
05700			BIG=-1.0E+20
05800			DO 360 J=1,NSP
05900			IF (X(J).EQ.0) X(J)=1.*10E-36
06000			SPT(J)=-10.*ALOG10(X(J))
06100	360		IF(BIG.LT.SPT(J))BIG=SPT(J)
06200			DO 377 J=1,NSP
06300	377		SPT(J)=SPT(J)-BIG
06400		RETURN
06500		END
     

00100		SUBROUTINE INVFL(X,NP,M,RO,ERRN)
00200		DIMENSION X(512)
00300		DIMENSION F(50),TF(50),A(50),R(50)
00400	C	CALCULATION OF M+1 LENGTH AUTOCORRELATION SEQUENCE
00500	C	FROM THE INVERSE FILTER INPUT SEQUENCE
00600		MP1=M+1
00700			DO 11 JJ=1,MP1
00800			J=JJ-1
00900			MMJ=NP-J
01000			SS=0.
01100				DO 10 I=1,MMJ
01200				IPJ=I+J
01300	10			SS=SS+X(I)*X(IPJ)
01400	11		R(JJ)=SS
01500		RO=R(1)
01600	C	FORTRAN IMPLEMENTATION OF LEVINSON'S METHOD FOR
01700	C	SOLVING THE AUTOCORRELATION MATRIX
01800		F(1)=1.
01900		ALPHA=R(1)
02000		BETA=R(2)
02100		A(1)=-R(2)/R(1)
02200		GAMMA=A(1)*R(2)
02300			DO 1 N=2,M
02400			NM1=N-1
02500			C=-BETA/ALPHA
02600			IF (N-2) 2,2,3
02700	3			DO 4 J=2,NM1
02800				NN=N-J+1
02900	4			TF(J)=F(J)+C*F(NN)
03000				DO 5 J=2,NM1
03100	5			F(J)=TF(J)
03200	2		F(N)=C*F(1)
03300			ALPHA=ALPHA+C*BETA
03400			BETA=0.
03500				DO 6 J=1,N
03600				NN=N-J+2
03700	6			BETA=BETA+F(J)*R(NN)
03800			Q=-(R(N+1)+GAMMA)/ALPHA
03900				DO 7 J=1,NM1
04000				NN=N-J+1
04100	7			A(J)=A(J)+Q*F(NN)
04200			A(N)=Q*F(1)
04300			GAMMA=0.
04400				DO 8 J=1,N
04500				NN=N-J+2
04600	8			GAMMA=GAMMA+A(J)*R(NN)
04700	1		CONTINUE
04800	C	CALCULATE NORMALIZED ERROR ERRN
04900	C	SM=0.
05000	C		DO 77 J=1,M
05100	C 77		SM=SM+A(J)*R(J+1)
05200	C	ERRN=1.+SM/R(1)
05300	C	PLACE COEFFICIENTS IN PROPER LOCATIONS OF THE
05400	C	REAL VECTOR X FOR FFT-ING TO GET INVERSE SPECTRUM
05500			DO 51 J=1,M
05600	51		X(J+1)=A(J)
05700		X(1)=1.
05800		MP1=MP1+1
05900			DO 123 J=MP1,512
06000	123		X(J)=0.
06100		RETURN
06200		END
     

00100		SUBROUTINE FFTP(X,Y,M,L)
00200		DIMENSION X(512),Y(512)
00300		N=2**M
00400		L2=2**L
00500		        DO 1 LO=1,M
00600		LMX=2**(M-LO)
00700		LMM=LMX
00800		LIX=2*LMX
00900		SCL=6.283185/LIX
01000	C     TEST FOR PRUNING
01100		IF (LO-M+L) 20,30,30
01200	20	LMM=L2
01300	30	        DO 1 LM=1,LMM
01400		ARG=(LM-1)*SCL
01500		C=COS(ARG)
01600		S=SIN(ARG)
01700		        DO 1 LI=LIX,N,LIX
01800		J1=LI-LIX+LM
01900		J2=J1+LMX
02000		T1=X(J1)-X(J2)
02100		T2=Y(J1)-Y(J2)
02200		X(J1)=X(J1)+X(J2)
02300		Y(J1)=Y(J1)+Y(J2)
02400		X(J2)=C*T1+S*T2
02500	1	Y(J2)=C*T2-S*T1
02600		RETURN
02700		END
     

00100		SUBROUTINE RBITS(X,Y,M)
00200		DIMENSION X(512),Y(512),L(9)
00300		EQUIVALENCE(L9,L(1)),(L8,L(2)),(L7,L(3)),(L6,L(4)),
00400	     1 (L5,L(5)),(L4,L(6)),(L3,L(7)),(L2,L(8)),(L1,L(9))
00500	C       PERFORM BIT REVERSAL
00600		    DO 70 J=1,9
00700		    L(J)=1
00800		    IF (J-M) 71,71,70
00900	71	    L(J)=2**(M+1-J)
01000	70	    CONTINUE
01100		JN=1
01200		    DO 60 J1=1,L1
01300		    DO 60 J2=J1,L2,L1
01400		    DO 60 J3=J2,L3,L2
01500		    DO 60 J4=J3,L4,L3
01600		    DO 60 J5=J4,L5,L4
01700		    DO 60 J6=J5,L6,L5
01800		    DO 60 J7=J6,L7,L6
01900		    DO 60 J8=J7,L8,L7
02000		    DO 60 JR=J8,L9,L8
02100		IF (JN-JR) 61,61,62
02200	61	R=X(JN)
02300		X(JN)=X(JR)
02400		X(JR)=R
02500		FI=Y(JN)
02600		Y(JN)=Y(JR)
02700		Y(JR)=FI
02800	62	JN=JN+1
02900	60	CONTINUE
03000		RETURN
03100		END
     

00100		SUBROUTINE USCRM(X,Y,M)
00200		DIMENSION X(512),Y(512)
00300		I=2**M
00400		LO2=I/2
00500		SA=3.141593/I
00600		X(1)=2.*(X(1)+Y(1))
00700		Y(1)=0.
00800		LLL=LO2+1
00900		X(LLL)=2.*X(LLL)
01000		Y(LLL)=-2.*Y(LLL)
01100			DO 33 K=2,LO2
01200			J=I-K+2
01300			ANG=SA*(K-1)
01400			C=COS(ANG)
01500			S=SIN(ANG)
01600			AA=X(K)+X(J)
01700			BB=Y(K)-Y(J)
01800			CC=X(K)-X(J)
01900			DD=Y(K)+Y(J)
02000			XR=C*DD-S*CC
02100			XI=S*DD+C*CC
02200			X(K)=AA+XR
02300			Y(K)=BB-XI
02400			X(J)=AA-XR
02500			Y(J)=-BB-XI
02600	33	CONTINUE
02700		RETURN
02800		END